% The code is based on the following paper:
%   Chen, T. Y., Lin, Y. L., and Tzeng, L. Y., forthcoming.
%   Estimating probability weighting functions through option pricing bounds. Review of Asset Pricing Studies.
%
% Copyright: Tzu-Ying Chen, Yo-Lan Lin, Larry Y. Tzeng
% Date: January 30, 2024

clear

%% Load Data
FileName = ['AllTradingDate19962021_Monthly_Pre1Day.mat'];
load(FileName, ...
     'Target_AllDate');             
clear FileName

%% Calculate At-the-Money Implied-Volatility 
ATMIV = nan(size(Target_AllDate, 1), 1);                                   
for d = 1:length(Target_AllDate)
    Target_Date = Target_AllDate(d, 1);
    Target_TTM = Target_AllDate(d, 3);

    % Load Data
    FileName = ['OP' num2str(fix(Target_Date / 10000)) '_' num2str(fix(rem(Target_Date, 10000) / 100)) '.txt'];
    Data = load(FileName);     
    clear FileName

    Index_ID = 1;
    Index_Date = 2;
    Index_TTM = 3;
    Index_CPFlag = 4;
    Index_K = 5;
    Index_S = 6;
    Index_OP_Bid = 7;
    Index_OP_Ask = 8;
    Index_OI = 9;
    Index_V = 10;
    Index_RF = Index_V + 1;                                                
    Index_DY = Index_V + 2;                                                
    Index_IS = Index_V + 3;                                                
    Index_IS_ADJ = Index_V + 4;                                            
    Index_S_ADJ = Index_V + 5;                                             
    Index_IV_Bid = Index_V + 6;                                            
    Index_IV_Ask = Index_V + 7;                                            

    Date_EXP_WeekDay = weekday(datenum(num2str(Data(:, Index_Date)), 'yyyymmdd') ...
                               + Data(:, Index_TTM));      
    Data(Date_EXP_WeekDay==7, Index_TTM) = Data(Date_EXP_WeekDay==7, Index_TTM) - 1;      
    clear Date_EXP_WeekDay            
    
    Index = find((Data(:, Index_Date)==Target_Date) & ...
                 (Data(:, Index_TTM)==Target_TTM));
    Data = Data(Index, :);                                                 
    Data(:, Index_TTM) = Data(:, Index_TTM) - 1;                           
    clear Index
    
    % Risk-Free Rate  
    FileName = ['RiskFreeRate19962021.txt'];
    Data_RF = load(FileName);
    clear FileName

    Data(:, Index_RF) = RF_TTM(Data_RF, Data(:, Index_Date), Data(:, Index_TTM));   
    clear Data_RF
    
    % Dividend Yield
    FileName = ['IndexDivYield19962021.txt'];
    Data_DY = load(FileName);
    clear FileName

    Data(:, Index_DY) = Data_DY(Data_DY(:, Index_Date)==Target_Date, end); 
    clear Data_DY
    
    % Implied Stock Price    
    K = Data(:, Index_K);
    OP = 0.5 * (Data(:, Index_OP_Bid) + Data(:, Index_OP_Ask));
    CPFlag = Data(:, Index_CPFlag);   
    S0 = Data(1, Index_S);
    TTM = Data(1, Index_TTM) / 365;                                                                          
    RF = Data(1, Index_RF);                                                    
    DY = Data(1, Index_DY);                                                   
     
    [AllK, ~, Index_AllK] = unique(K);                           
    Pair_OP_AllK = accumarray([Index_AllK CPFlag], ...
                              OP, ...
                              [length(AllK) 2], ...
                              @mean, ...
                              nan);                           
    Index_Drop = find(sum(isnan(Pair_OP_AllK), 2) > 0);
    if length(Index_Drop) > 0
        AllK(Index_Drop) = [];                                             
        Pair_OP_AllK(Index_Drop, :) = [];                                  
    else
    end
    clear Index_Drop   
    
    [~, Index_ATM] = min(abs(AllK / S0 - 1));
    K_ATM = AllK(Index_ATM);
    OP_C_ATM = Pair_OP_AllK(Index_ATM, 1);
    OP_P_ATM = Pair_OP_AllK(Index_ATM, 2);
    Data(:, Index_IS) = (OP_C_ATM - OP_P_ATM + K_ATM * exp(- RF * TTM)) / exp(- DY * TTM);    
    clear Index_AllK AllK Pair_OP_AllK Index_ATM K_ATM OP_C_ATM OP_P_ATM

    Data(:, Index_IS_ADJ) = exp(- DY * TTM) * Data(:, Index_IS);  
    Data(:, Index_S_ADJ) = exp(- DY * TTM) * Data(:, Index_S); 
    clear K OP CPFlag S0 TTM RF DY

    % Data Filtering 
    K = Data(:, Index_K);                                                                                          
    S0 = Data(:, Index_IS);  
    Index_C = find((Data(:, Index_CPFlag)==1) & ...
                   ((K ./ S0) >= 1) & ...
                   ((K ./ S0) <= 1.05));    
    Index_P = find((Data(:, Index_CPFlag)==2) & ...
                   ((K ./ S0) >= 0.9) & ...
                   ((K ./ S0) < 1));                
    Index = union(Index_C, Index_P);
    [~, Index_Min] = min(abs((K(Index) ./ S0(Index)) - 1));
    Data = Data(Index(Index_Min), :);                                      
    clear Index Index_C Index_P Index_Min K S0
        
    % Option-Implied Volatility 
    Type = cell(size(Data, 1), 1);                                         
    Type(Data(:, Index_CPFlag)==1) = {'Call'};   
    Type(Data(:, Index_CPFlag)==2) = {'Put'};
    K = Data(:, Index_K);
    OP = 0.5 * (Data(:, Index_OP_Bid) + Data(:, Index_OP_Ask)); 
    S0_ADJ = Data(:, Index_IS_ADJ);   
    TTM = Data(:, Index_TTM) / 365;                                        
    RF = Data(:, Index_RF);                                                    
    ATMIV(d) = blsimpv(S0_ADJ, K, RF, TTM, OP, [], [], [], Type);             
    clear Type K OP S0_ADJ TTM RF       

    clear Target_Date Target_TTM Data
    clear Index_ID Index_Date Index_TTM Index_CPFlag Index_K Index_S Index_F Index_OP_Bid Index_OP_Ask Index_OI Index_V
    clear Index_RF Index_DY Index_IS Index_IS_ADJ Index_S_ADJ Index_IV_Bid Index_IV_Ask
end
clear d 

%% Output
FileName = ['ATMIV19962021_Monthly_Pre1Day'];
save(FileName);
clear FileName 